Method of Segmenting Anatomic Entities in Digital Medical Images

ABSTRACT

For each of a number of landmarks in an image an initial position of the landmark is defined. Next a neighborhood around the initial position, comprising a number of candidate locations of the landmark is sampled and a cost is associated with each of the candidate locations. A cost function expressing a weighted sum of overall gray level cost and overall shape cost for all candidate locations is optimized. A segmented anatomic entity is defined as a path through a selected combination of candidate locations for which combination the cost function is optimized.

RELATED APPLICATIONS

This application claims priority to European Application Nos. EP05107903.6, filed Aug. 30, 2005, and EP 05107907.7, filed Aug. 30, 2005,and also claims the benefit of U.S. Provisional Application No.60/715,878, filed on Sep. 9, 2005, all three of which are incorporatedherein by reference in their entirety.

FIELD OF THE INVENTION

The present invention relates to a method of segmenting an entity in adigital image, more specifically an anatomic entity in a digital medicalimage. The segmentation process according to the present invention isbased on the use of geometric and photometric models of an image.

The method of the present invention can inter alia be used as a processof performing geometric measurements on digital images, morespecifically on digital medical images.

BACKGROUND OF THE INVENTION

In radiological practice, geometric measurements are frequently used toaid diagnosis of abnormalities. In order to perform these measurements,key user points must be placed in an image, for example in an imagedisplayed on a display device, on their corresponding anatomicallandmark position. Measurements such as the distance between two points,or the angulation between lines are based on the position of these keyuser points. Furthermore, the geometry as a whole may be assessed fornormality or abnormality, involving an analysis of the complete shape.Hence there is a need to automate and objectify the extraction ofquantitative information that is embedded in a radiological image.

An example of such a frequently performed measurement is the computationof the cardiothoracic ratio (CTR) in thorax RX images (FIG. 6). Thisratio is defined as the ratio of the transverse diameter of the heart,at the level of the apex, to the internal diameter of the thorax (ID),i.e. CTR=(MLD+MRD)/ID.

The transverse diameter of the heart is composed of the maximumtransverse diameter on the left side of the heart (MLD) and the maximumtransverse diameter on the right side of the heart (MRD). Clearly, thisdefinition entails that the radiologist searches along the inner borderof the left and right ribcage boundary to locate the point pair neededto compute the internal diameter ID. This point pair must lie at thegreatest internal diameter of the thorax. Likewise, the left and rightheart shadow border must be searched to locate the points needed tocompute the sum of MLD and MRD. More specifically, these points aresituated most distant with respect to the midline of the spine. Theprocess of border search requires that the radiologist is performinganatomy segmentation and locating the points (a total of four in thisexample) on the segmented anatomy. The segmentation step, in the case ofCTR computation, amounts to delineating the lung fields.

Many other measurements in digital images follow a similar approachinvolving the segmentation (224,418) of the anatomic organ or entityonto which segmented geometry the characteristic points (510) andmeasurement objects (512) are determined and finally measurements (514)are performed (FIG. 5).

Referring to the example of cardiothoracic index calculation, toautomatically position the required points, a method is needed toautomatically segment the lung fields on a chest radiographic image.

The segmentation problem can be approached in several ways, depending onthe application. Segmentation strategies evolved from low-levelstrategies in the early years of computer vision to the more recentmodel-based strategies.

Low-level methods rely on local image operators separating pixels withdifferent photometric characteristics and grouping of pixels withsimilar local photometric characteristics. Examples of both classes areedge detection and region growing. Despite the poor performance of theselow-level approaches, they are very popular in most commercial imageanalysis tools. The main reasons are that they are simple to understandand to implement. For complex image data however, such as present inmedical images and exemplified by the content of a thorax image asdescribed above, their usefulness is limited.

More successful methods incorporate a priori knowledge about the shapeto be segmented and about the photometric or gray-level appearance ofthe object in the image. These methods, referred to as model-basedmethods are often based on template matching. A template is matched forinstance by correlation or with generalized Hough transform techniques.Unfortunately, the template matching is likely to fail in case ofmedical images. This is due to the large variability in shape andgray-level appearance that the anatomic object may exhibit.

Methods based on active contours, introduced by Kass et. al. (M. Kass,A. Witkin, and D. Terzopoulos, Snakes: active contour models, Int. J.Computer Vision, 1(4):321-331, 1988) and level sets (J. A. Sethian,Level set methods and fast marching methods, Cambridge Univ. Press,Cambridge, U.K. 1999) are able to cope with a larger shape variability,but are still unsuited for many medical segmentation tasks becauselittle a priori knowledge about the object to be segmented can beincorporated. Handcrafted parametric models overcome this problem, butare limited to a single application.

In view of these shortcomings, it is obvious that there is need for ageneric segmentation scheme that can be trained with examples in orderto acquire knowledge about the shape of the object to be segmented andthe gray-level appearance of the object in the image. Active shapemodels (ASMs), introduced by Cootes and Taylor (T. F. Cootes, C. J.Taylor, D. Cooper, J. Graham, Active Shape Models—their training andapplications, Computer Vision and Image Understanding, 61(1):38-59,1995) satisfy this definition of segmentation schemes. The shape modelis given by the principal components of vectors of landmark points. Thegray-level appearance model describes the statistics of the normalizedfirst derivative of profiles centered at each landmark that runperpendicular to the object contour. The location of a landmark in a newimage is found by minimizing the Mahalanobis distance between the firstderivative profile and the distribution of the profile. This algorithmstarts from an initial estimate and performs a fitting procedure, whichis an alternation of landmark displacements and shape model fitting.Similar approaches have been devised all employing a three-stepprocedure. First, they all use a shape model that ensures that plausibleresults are generated. Secondly, they use a gray-level appearance modelto place the object at a location where the gray-level pattern aroundthe border or within the object is similar to what is expected from thetraining examples. Finally, the algorithm fits the model by minimizingsome cost function.

SUMMARY OF THE INVENTION

The approach presented by active shape models is still faced withseveral limitations.

A first limitation is the need for an initial estimate, known as theinitial positioning problem. In some cases, an extensive search for aproper initial contour is needed in this framework.

A second limitation lies in the alternate use of shape model andgray-level appearance model. Firstly, the segmentation estimate isupdated using the gray-level appearance model. Secondly, the shape modelis fitted to the updated contour. Unfortunately the shape model ismisled if the gray-level appearance model wrongly locates a landmark.

Some experiments with active shape model segmentation schemes exhibitedanother problem. If the gray-level appearance model for a specificlandmark is used to find a better location for that landmark, itrequires that the true landmark location lies inside the region thealgorithm is exploring. In some cases the algorithm may search for alandmark in the region of another landmark. This means that the wronggray-level appearance model is used at the wrong region, resulting in apoor segmentation.

It is an object of the present invention to provide a method ofsegmenting an entity in a digital image, more specifically an anatomicentity in a digital medical image, that overcomes the problems of theprior art.

The above-mentioned aspects are realized by a method as set out in claim1.

Specific features for preferred embodiments of the invention are set outin the dependent claims.

Another aspect of this invention relates to a method of measuring ananatomic entity in a digital medical image as set out in the appendingclaims.

In the following the terms gray value model, gray value appearancemodel, gray level model, gray level appearance model and photometricmodel are used as synonyms.

Likewise the terms shape model and geometric model are used as synonyms.

The terms feature image and feature are also used as synonyms.

The embodiments of the methods of the present invention are generallyimplemented in the form of a computer program product adapted to carryout the method steps of the present invention when run on a computer.The computer program product is commonly stored in a computer readablecarrier medium such as a DVD or a CD-ROM. Alternatively the computerprogram product takes the form of an electric signal and can becommunicated to a user through electronic communication.

In EP A 1349098, which is incorporated herein by this reference, amethod is disclosed to automate the measurements in digitally acquiredmedical images by grouping measurement objects and entities into acomputerized measurement scheme including a bi-directionally linkedexternal graphical model and an internal informatics model. In ameasurement session according to EP A 1349098, a measurement scheme isretrieved from the computer and activated. Measurements are subsequentlyperformed on the displayed image under guidance of the activatedmeasurement scheme.

In this computerized method, a multitude of geometric objects are mappedin the digital image onto which other measurement objects and finallymeasurement entities (such as distances and angles) are based. The basicgeometric objects are typically key user points, which define othergeometric objects onto which geometric measurements are based. The aboveidentified patent application does not disclose however how the mappingmay be effectuated automatically without the need of user positioningand manipulation of key measurement points. The method according to thepresent invention can be used to effectuate this mapping.

According to the present invention, a segmentation method is providedwhich comprises the construction of models of anatomic entities such asan organ or a structure in a medical image, which models can be trainedfrom examples so as to obtain knowledge of the gray value appearance andknowledge of the shape of the anatomic entity to be segmented.

The gray value appearance model comprises the probability distributionof the gray value and a number of multiscale derivatives of gray valuesat a set of landmarks distributed over the anatomic outline.

The shape model comprises the probability distribution of eachconnection vector between successive landmarks along the outline of ananatomic entity.

A discrete point-based object representation is introduced to describethe anatomic outline in a medical image.

Two strategies and associated systems for generating models andapplication of the models to segment actual images are disclosed.

These strategies, although decomposing in building blocks with equalgoals, construct and employ different photometric and geometric modelknowledge.

First, a new gray level appearance model is constructed from trainingimages, encoding photometric knowledge at landmark positions. This stepexploits intensity correlation in neighborhoods sampled around eachlandmark.

Secondly, a shape model is constructed, encoding geometric knowledgebetween landmarks. This step exploits spatial correlation betweenlandmarks of the segmented objects.

In a segmentation step, photometric and geometric knowledge are jointlyused to segment one or more anatomic contours on a new image.

The resulting segmentations may be used to derive the position ofdesired measurement points in the course of performing geometricmeasurements on an image.

Measurement points can be defined in either a direct or an indirect way.

In the direct way, the positions of the landmarks are used as individualmeasurement points.

In the indirect way, the landmarks resulting from the segmentationprocess are interconnected by a curve and characteristic points on thiscurve are defined as measurement points.

Gray Level Appearance Models

The first step in each modeling system constrains the geometric positionof the landmark to lie on a linear path perpendicular to the presentcontour or a sparse grid centered on the landmark. The photometricpattern at a landmark point in an image is captured by a profile vectorsampled in feature images based on image derivates extracted at a numberof scales. The profile sampling may be along a linear path in the imagesor on a circular path around the landmark point.

Feature Images

When the photometric pattern is computed for all image points in theimage, feature images are obtained, the size of which is equal to thatof the original image. Neighborhood operations are typically employed tocompute the feature images. The values of the features in the featureimage for a specific point can be arranged in a feature vector. Thecomputation of feature images typically comprise two steps: (a) a linearfiltering step using a bank of convolution kernels and (b) apost-processing step involving the sub-steps of a non-linear pointoperation combined with the computation of local image statistics. Instep (b) one of the sub-steps may be performed alone.

Different filters and different types of post-processing lead todifferent types of feature images. For example, Laws constructed 25two-dimensional convolution kernels as the outer product of pairs ofvectors resembling ideal feature profiles such as Level, Edge, Spot,Wave and Ripple. The post-processing step comprises non-linear windowingin which the absolute filtered values are averaged in larger windows toobtain image texture energy measures. Unser used filter banksconstructed as the outer product of vectors corresponding to well-knowntransforms like discrete sine (DST), discrete cosine (DCT), discreteeven sine (DEST), discrete real even Fourier (DREFT), discrete real oddFourier (DROFT) transforms, and for post-processing the computation ofthe channel variances. Finally, Gabor filter banks are constructed usingsymmetric or anti-symmetric Gabor filters tuned to combinations ofdifferent spatial frequencies and different orientations, such that alloutputs have equal frequency bandwidth. The non-linear point operationstep may include thresholding in this case.

Other image characteristics may be measured using feature images basedon multi-resolution wavelets filters, rank-value filters or adaptivefilters.

Feature images based on locally orderless images (LOI's) The imagestructure around an anatomical landmark in a medical image exhibits atypical gray-level pattern. In the present invention, this structure iscaptured in a number N of mathematical features as outlined in thesequel. The concept of using feature images to build a gray-levelappearance model was introduced earlier in B. Van Ginneken et al.,Active shape model segmentation with optimal features, IEEE Trans. onMedical Imaging, 21(8):924-933, 2002, based on so-called locallyorderless images (LOI's) presented by J. J. Koenderink and A. J.Vandoorn, The structure of locally orderless images, Int. J. of ComputerVision, 31(2):159-168, 1999.

A Taylor expansion can be used to approximate a gray-level function ƒaround a point of interest, i.e. a landmark, at position x₀, by apolynomial of order K. The coefficients of each term are given by thederivatives ƒ^((i)) at x₀:${f(x)} \approx {\sum\limits_{i = 0}^{K}{\frac{1}{i!}{f^{(n)}\left( x_{0} \right)}\left( {x - x_{0}} \right)^{i}}}$

All information in this function approximation is contained in thederivatives ƒ^((i)). Similarly, the image can be described by an imagefilter bank containing the derivatives of the original image(L,L_(x),L_(y),L_(xx),L_(yy),L_(xy), . . . ). Each derivative image mayfurther be computed at a number of blurring or inner scales σ. Thefeature images are finally obtained by computing the post-processingstep as local image statistics of each derivate image. The local imagestatistics comprise a number of moments of the local histogram in awindow with width α around each location x₀. In a specific embodimentthe first and second histogram moments, all derivatives up tosecond-order (L,L_(x),L_(y),L_(xx),L_(yy),L_(xy)) and 5 inner scales(σ=0.5,1,2,4,8 pixels) are computed, amounting to a total of N=2×6×5=60feature images. The window scale to compute the histogram is inaccordance to the inner scale, i.e. α=2σ.

The limitation of using the derivatives in the row and column directly(i.e. x resp. y) in the feature vector is the lack of invariance withrespect to translation and rotation (i.e. Euclidean Transformations).Unless the images are anti-rotated into a standard orientation, theseoperators can only be trained on examples that have the sameorientation. Cartesian differential invariants (CDI) describe thedifferential structure of an image independently of the chosen Cartesiancoordinate system; hence the anti-rotation step can be omitted. Theconstruction of CDI's from Gaussian differential operators was describedin L. M. J. Florack et al., Scale and the differential structure ofimages, Image and Vision Computing, Vol. 10 (6):376-388, 1992. Thefollowing are CDI's using derivatives up to second order:(L,L_(xx)+L_(yy),L² _(x)+L² _(y),L² _(x)L_(xx)+2L_(xy)L_(x)L_(y)+L²_(y)L_(yy),L² _(xx)+2L² _(xy)+L² _(yy)). Again, at each pixel location,feature images are constructed using these operators at number of innerscales σ and computing the first two moments of the locally orderlesshistograms in windows of extent α, with α linked to the inner scale,e.g. α=2σ. These operators are independent to Euclidean transformationsof image geometry, but still depend on the choice of the scale parameterσ.

Profile and Feature Similarity Measurement

The gray level appearance model is built on the basis of a mean profileand variance-covariance matrix, and exploits intensity correlationbetween points of the profile, as explained next. The intensity at aprofile is characterized by the intensity feature vector as outlinedabove.

In order to compare a profile with another profile, it is necessary toknow how each point in the profile relates to any other point, and tomeasure the extent by which the value of a certain point of the profilevaries from the mean with respect to any of the other points. Thisstatistical relationship is captured by the covariance measure, which isalways computed between 2 points (or the features at those points). Ifthe covariance is measured between the feature value of a point anditself, the variance is obtained. When the profile includes k points ateither side of the contour, for a total of 2k+1 points, (2k+1)k numberof relationships exists. The covariance between a pair of points fromthe profile tells the way how a first value of the pair varies withrespect to the second. When the covariance is large and positive, itindicates that when the first point's feature value increases, so doesthe second in a proportional way. When the covariance is large andnegative, it indicates that when the first point's feature valueincreases, the second decreases in a proportional way. When thecovariance is zero, it means that there is no systematic coupling in thefeature values, i.e. they are independent of each other. Allvariance-covariance pairs can be arranged in a covariance matrix S,which is symmetrical about the main diagonal. The covariance matrix thuscaptures all structural dependencies between feature value pairs in theprofile. In practice, the covariance matrix is built from learningsamples, in this case a collection of profiles at a number of landmarksand a number of intensity features in a set of training images.

The covariance relates to how distance is calculated in higherdimensional space, i.e. the space of the 2k+1 points in a given profile,to compute the distance between a current profile and a model profile.When the profiles would only comprise a single point with value x, theprofiles could be compared with the model point value x for similarityjust by subtracting the values of the current profile point and modelprofile point and taking the absolute value, i.e. d=|x− x|. When theprofiles would comprise two points, the distance between the 2-vector xand the model 2-vector x could be computed as the Euclidean distance inthe plane, i.e.d=√{square root over ((x− x)′(x− x )=√{square root over (|x₁− x ₁|²+|x₂−x ₂|².Euclidean distance weighs all directions equally. However, when thevariables x₁ and x₂ have unequal standard deviations σ₁ resp. σ₂, thecontours of iso-distance lines are ellipses with major axes parallel tothe coordinate axes. The differences may be weighted with the inverse ofthe standard deviation, resulting in the weighted Euclidean distance$d_{w} = {\sqrt{\left( {x - \overset{\_}{x}} \right)^{\prime}{W\left( {x - \overset{\_}{x}} \right)}} = {\sqrt{{\frac{x_{1} - {\overset{\_}{x}}_{1}}{\sigma_{1}}}^{2} + {\frac{x_{2} - {\overset{\_}{x}}_{2}}{\sigma_{2}}}^{2}}.}}$The matrix W contains the inverse of the variances of the variables,i.e. $W = {\begin{bmatrix}{1/\sigma_{1}^{2}} \\{1/\sigma_{2}^{2}}\end{bmatrix}.}$Intensity levels of pixels close to each other as in the case of linearor circular profiles will be heavily correlated when they belong to anidentical anatomical region with smoothly varying intensity levels. Whenpixel pairs are considered belonging to different tissues, thecorrelation may be inversely proportional. In either case, imaging noiseand unimportant anatomical variation introduces a non-zero variance forthe intensity level of the pixel itself.

Hence, if the variables are also correlated and have differentvariances, the inverse of the covariance matrix must be inserted,yielding the so-called Mahalanobis distanced _(m)=√{square root over ((x− x )′S ⁻¹(x− x ).

The Mahalanobis distance weights the distance of a multi-dimensionaldata point x from its mean x such that observations that are on the samemultivariate normal density contour will have the same distance d_(m).These Mahalanobis iso-distance contours take the form of ellipses intwo-dimensional space, and ellipsoids in three-dimensional space.Obviously, the number of points in the profile 2k+1 is normally muchhigher than three, in which case the iso-distance loci take the form ofhyper-ellipsoids in multi-dimensional space.

The gray level appearance model is obtained from a training image set,and includes for each landmark in a mean profile and a covariance matrixfor each of the features. The total of stored models is thusproportional to the number of landmarks and the number of features.

Shape Models

The positional correlation between landmarks along the contour is takeninto account in a shape modelling step. Because segmentation of anatomicstructures with relatively deterministic shape is aimed at, knowledgeregarding the shape must preferably encode the deterministic component,along with relevant anatomical variation, and further preferably ignoreirrelevant shape details.

Two embodiments are outlined to take into account positionalcorrelation. The first is based on a global approach involving the useof all landmark positions simultaneously; the second employs a localapproach, modelling the spatial disposition of successive landmarks.

PCA Analysis and Fitting

In the first embodiment, a variance-covariance matrix of landmarkcoordinate positions is constructed from manually placed landmarks in aset of aligned training images. Each diagonal entry in thevariance-covariance matrix represents the variance of the coordinatevalues of the landmarks. The non-diagonal entries represent theco-variance of pairs of landmark coordinate values, i.e. it expresseshow the coordinates of two landmarks vary with respect to each other.The variation may systematically be in equal ‘positive directions’ withrespect to their means, yielding high positive covariance, orsystematically in equal ‘negative directions’ with respect to theirrespective means, yielding high negative covariance, or randomly inopposite directions with respect to their means, yielding zerocovariance. An eigenvector/eigenvalue analysis, known in the prior artas Principal Components Analysis (PCA), will detect any globalsystematic variation modes of the landmarks with respect to the meanlandmark positions. PCA aims to reduce the dimensionality of a data setby only keeping the components of coordinate set with largest variation;it therefore statistically decorrelates the original coordinatepositions. By projecting onto these highly varying subspaces, therelevant statistics can be approximated by a smaller dimensional systemthat retains important shape features. A sufficient number of theseprincipal modes are retained and used as the shape model, ignoring theother modes attributed to noisy variations.

Landmark Connection Vector Model

The second embodiment aims at modeling the disposition of successivelandmarks at the local level. However, this restriction does not implythat the global shape is not captured. The shape outline can becompletely reconstructed given the start point coordinates and theconnection vectors, by iteratively adding the connection vector to thecurrent point, starting from the start point. This property is exploitedin shape description at the pixel level by chain codes for example. Theconnection vector is a discrete approximation of the local tangentvector to the curve. When taking into account the direction vectorbetween successive landmarks, i.e. first order geometric derivatives, itallows describing the shape of the object irrespective of translationsof the object. By taking into account the curvature, i.e. second ordergeometric differentiation, between successive landmarks, for exampleexpressed by the angle subtended by two successive connection vectors,the shape of the object can be described irrespective of translationsand rotations of the shape (so-called rigid body transformations). Inthat case the curve can be completely reconstructed when the startposition, and the first geometric derivative (i.e. start connectionvector) is known, by iteratively adding the difference vector ofsuccessive tangent vectors. Curvature based on a chain code at the pixellevel is usually highly noise dependent for the scale is at a too finelevel. This drawback is coped with in the present disclosure because thelandmarks are sufficiently spaced apart so as to enable computation ofnoise-independent yet accurate first or second order derivativeinformation. Because the magnitude of tangent or curvature vectordepends on scale, normalizing the size of each curve of the training setto unity for example, further allows describing the curve up to anEuclidean similarity transformation. The mean and covariance statisticsof displacement vectors (i.e. first derivative or tangent) or vectordifferences of displacement vectors (i.e. second derivative orcurvature) between manually placed landmarks in a set of alignedtraining images constitutes the shape model of this embodiment.

Segmentation Systems

Segmentation systems can also be organised along two tracks each havingbuilding blocks with related objectives.

A first building block aims at confining the admissible positions ofinstantiated model landmark positions to a number of most probablelocations in a target image.

In the case landmarks are represented by linear profiles in featureimages, a best point voting approach can be adopted that ranks a set ofcandidate positions on a perpendicular line through an instantiatedlandmark.

For a circular profile case, a rectangular search grid can be consideredat each instantiated landmark and grid crossing positions are rankedaccording to the combined Mahalanobis distance of circular profilesthrough all feature images to the model circular profiles. Second andthird building blocks jointly constitute an optimization step. A costmatrix may be constructed, using the cost measures associated with theretained admissible position for each landmark. The cost comprises aterm associated with agreement in gray level appearance. Agreement withthe shape model may be verified at the global level using weighted PCAcontour fitting. Alternatively, the shape agreement may be verified byincorporating a shape cost associated with the candidate position of thelandmark. A dynamic programming algorithm can be used to construct aminimum cost path through the cost matrix, yielding the finalsegmentation as the path with lowest overall cost.

Finally, characteristic points on the final segmentation may indirectlydefine the geometrical position of the key measurement points, and alldepending measurement objects and measurement quantities may begenerated on the basis of these key measurement points.

Alternatively, the landmark points used in the segmentation may be usedto directly define the key measurement points and a search grid for eachmeasurement point may be provided to locate the admissible positions ofthe measurement point.

Further specific embodiments and associated advantages of the presentinvention will become apparent from the following description anddrawings.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a flow chart illustrating a first modeling system (Modelingsystem I) according to the present invention,

FIG. 2 is a flow chart illustrating a first segmentation system(Segmentation System I) according to the present invention,

FIG. 3 is a flow chart illustrating a second modeling system (ModelingSystem II) according to the present invention,

FIG. 4 is a flow chart illustrating a second segmentation system(Segmentation System II) according to the present invention,

FIG. 5 is a flow chart illustrating measurement computation on the basisof the characteristic points on the anatomic objects segmented accordingto the present invention,

FIG. 6 is a measurement scheme for measurement of the cardiothoracicratio from a chest radiograph,

FIG. 7(a) illustrates manual delineation of the left lung in a chestradiograph, in FIG. 7(b) the left lung is represented as a set oflandmarks {p_(i)}^(n) _(i−1),in this example n=14.

FIG. 8 illustrates a discrete representation of the lungs fields by aset of n landmarks p₁=(x₁,y₁), . . . , p_(n)=(x_(n),y_(n)). The imagesizes are normalized between 0 . . . 1. The continuous curve isinterpolated through the discrete landmarks.

FIG. 9 illustrates a linear and a circular gray-level profile as asampling of gray values at a number of points on a line resp. a circlewith radius r_(c) around the landmark,

FIG. 10 shows experimentally estimated distributions for some connectionvectors on the left lung shape having 14 landmarks: (a) v_(l), (b) v₂,(c) V₃.

FIG. 11 shows first moment feature images of a chest radiograph,corresponding to the derivatives up to second-order(L,L_(x),L_(y),L_(xx),L_(yy),L_(xy)) and the scales σ equal to 0.5, 2and 8 pixels,

FIG. 12 shows second moment feature images of a chest radiograph,corresponding to the derivatives up to second-order(L,L_(x),L_(y),L_(xx),L_(yy),L_(xy)) and the scales σ equal to 0.5, 2and 8 pixels,

FIG. 13 shows examples of application of the present invention to thesegmentation of lung fields on a chest radiograph: (a) manualsegmentation of left and right lung field; (b) automatic segmentation ofleft and right lung field.

FIG. 14 illustrates the computation of the Cardiothoracic Ratio using(a) the manual lung field segmentations resulting in a CTR_(man)=0.47;(b) using the automatically derived lung field segmentations accordingto the present invention resulting in a CTR_(automatic)=0.45.

DETAILED DESCRIPTION OF THE INVENTION

The present invention will be explained in detail with regard to aspecific application, namely the segmentation of the lung field in amedical image.

Object Representation

In the specific embodiments of the method of the present invention thatare described below, an anatomical object in an image is representedmathematically as a fixed number of discrete, labeled points lying onthe contour that encloses the object, i.e. p₁=(x₁,y₁), . . . ,p_(n)=(x_(n),y_(n)).

The contour {p_(i)}^(n) _(i−1) runs from p₁ to p_(n) and back to p₁.Hence, the object may be captured by a discrete shape vector x=(x₁,y₁, .. . x_(n),y_(n))^(T). The coordinate system is chosen such that allpoints inside the image area lie in the domain [0,1]×[0,1] (FIG. 7).

The segmentation scheme described below needs a number of trainingimages, in which the shape vectors x are manually determined. Once thealgorithm is trained on the data set, it is able to produce the shapevector in a new image.

The continuous contour of an anatomical object can be reconstructed bycurve interpolation from the points {p_(i)}^(n) _(i=1). A first orderinterpolation links successive points with a line segment, resulting ina polygonal approximation. Higher order interpolations may be used toachieve a smoother outline. Irrespective of the order, the continuousrepresentation can be used to derive spatial properties of the curve,such as the tangent direction along the curve, or the normal directionat a certain point on the curve.

Modeling System I (FIG. 1)

Gray Level Appearance Model

The gray level appearance model (118) (also called photometric model)captures the spatial intensity distribution in the neighborhood ofso-called landmarks (also called landmark points). The spatial componentof the gray level appearance model has been devised along the use ofgray-level profiles. As will be apparent from their construction, otherschemes can be employed.

The number n of landmarks (110) that is for example used to accuratelydescribe a lung field outline may be as low as 14, e.g. landmark 1denotes the uppermost point of each lung field, landmark 7 is positionedat the cross section of either the left or right heart shadow with theirrespective hemi-diaphragm, landmark 9 denoting the costophrenic angle.The other landmarks can be equidistantly positioned between these majorlandmarks. Obviously, a higher number of landmarks may be used.

Linear Gray-Level Profiles (FIG. 9)

At each landmark, the gray level appearance model describing the typicalimage structure around the landmark is sampled perpendicular to thecontour. On either side of the landmark, k pixels are sampled using acertain step size, which gives profile vectors (112) of length 2k+1.This sampling scheme has the advantage that it is able to model linearintensity gradients at the landmark. The profiles may be normalized sothat the sum of absolute values of elements in the profile is 1. Thedirection of sampling can be computed as the mean normal directionbetween the normals on the line segments resulting from connecting thelandmark with its previous or its next neighbour, or as the normal onthe line segment connecting the previous and next landmark.

Gray-Level Feature Images

The image structure around a landmark exhibits is captured in a number Nof mathematical features as explained higher. The feature images (114)are obtained by computing a number of moments of the local histogram ofderivative images in a window with width α around each location x₀. In apreferred embodiment the first and second histogram moments, allderivatives up to second-order (L,L_(x),L_(y),L_(xx),L_(yy),L_(xy)) and5 inner scales (σ=0.5,1,2,4,8 pixels) are computed, amounting to a totalof N=2×6×5=60 feature images. The window scale to compute the histogramis in accordance to the inner scale, i.e. α=2σ. The set of featureimages (114) can be regarded as a gray-level decomposition of thegray-level function in the neighbourhood of a landmark.

An example of the computed feature images in a thorax image is given inFIG. 11 (first histogram moments of original image and 5 derivativesimages at three scales) and FIG. 12 (second histogram moments oforiginal image and 5 derivatives images at three scales). Thorax imagesnormally appear in the upright position in the image; anatomic objectspertaining to other radiographic examinations may exhibitnon-standardized position and rotation, in which case the constructionof feature images may employ Cartesian Differential Invariants.

Gray-Level Appearance Model

All images of a training set are used to build a statistical model foreach feature and each landmark. Denoting the normalized profile sampledat a landmark i (e.g. 14 in total for the 2D case of a thorax lungfield) in feature j (e.g. 60 in total for 2D images) as a vectorg_(i,j), the probability distribution of g_(i,j) is estimated from thetraining images, resulting in a mean g _(i,j), and a covariance S_(i,j)(116). For a length of a profile vector of 2k+1, the size of thiscovariance matrix (116) is (2k+1)×(2k+1). The gray level appearancemodel for a certain landmark i having the mean profiles g _(i,j),j=1 . .. N and covariance matrices S_(i,j), j=1 . . . N. The gray-levelappearance model for the total lung field having the collection of allmean profiles g _(i,j),i=1 . . . n,j=1 . . . N and covariance matricesS_(i,j),i=1 . . . n,j=1 . . . N.

Gray-Level Cost

The Mahalanobis distance between a new profile g_(i,j) sampled at apoint p in a new feature image j for a certain landmark i is computed ash_(i,j)(p)=(g_(i,j)(p)− g _(i,j))^(T)S⁻¹ _(i,j)(g_(i,j)(p)− g _(i,j))

A smaller Mahalanobis distance means a larger probability that theprofile g_(i,j)(p) originates from the Gaussian distribution with mean g_(i,j) and covariance S_(i,j). Hence, the Mahalanobis distance may beused as a gray-level cost measure, denoted as h_(i,j)(p). This costmeasure is a function of the location p in a given image. The locationthat most likely is the true location for the landmark i, in accordanceto feature j, is the point p that minimizes h_(i,j)(p). This gray-levelcost can thus be defined for each feature j.

Shape Model

The curve outline, in the case of lung fields: one for each of the lungfields, is described by n landmark points (110) (FIG. 8). These landmarkpoints are manually determined in a set of s training images, yielding asequence of points (x₁,y₁) . . . (x_(n),y_(n)). These coordinate tuplesare subsequently positioned in a vector x=(x₁,y₁, . . . ,x_(n),y_(n))^(T), representing the curve. Next Principal ComponentsAnalysis (PCA) (120) is applied to the shape vectors x of the trainingimages. The PCA projects the curve in a lower dimensional space,covering the most important modes of geometric variation of the lungfields. Each curve x ∈

can be approximated by b ∈

with t<<2n:x≈ x+Φbwith x the mean shape (122) and Φ ∈

the eigenvectors of the covariance matrix of x corresponding to the tlargest eigenvalues. Each eigenvalue determines the variance of theshape for that shape mode (124). This curve approximation, representedby the vector b, constitutes the shape model (126) and fits the curve xinto the shape model. The eigenshapes can be regarded as a zero-order(positional) geometric decomposition of the shape represented by the setof landmarks. The role of the shape model is to constrain thedeformation of the object between the bounds of variation imposed by thetraining set, for example the three standard deviations with respect tothe mean shape.Segmentation System I (FIG. 2)

The segmentation algorithm will initially position a curve, representinga lung field outline, in the image, and finds a segmentation solution(224) by iteratively updating the positions of the curve's constituentlandmarks until an optimization criterion reaches an optimal value. Eachiteration includes the following steps.

Step 1. Best Point Voting

Assuming a current position of a particular landmark, a better locationfor this landmark will be searched on a line perpendicular to thecontour. On each side of the current position of the landmark, n_(s)positions are evaluated. One of the 2n_(s)+1 points will be voted as thebest next point. Each feature of the set of N features contributes onevote to one of the 2n_(s)+1 points. That point of the 2n_(s)+1 points isselected that has minimum cost according to the Mahalanobis distancecriterion (214). For each of the 2n_(s)+1 points, the profile having2k+1 points sampled in feature image j is put in the equation ofh_(i,j)(p) with the proper mean g _(i,j) and covariance S_(i,j) filledin. The point with the lowest Mahalanobis distance is voted by thatfeature (216). Each feature chooses one point, resulting in N votesdivided over 2n_(s)+1 points as sets N₁ , . . . , N_(2n) ^(s) ₊₁ with${\sum\limits_{k = 1}^{{2n_{s}} + 1}N_{k}} = {N.}$Obviously the number N of features must be substantially larger than thenumber 2n_(s)+1 of selectable points for otherwise there are too fewvotes that may be assigned to a point of the set of selectable points.

The number of votes N_(i,j) that each point of the set of considered2n_(s)+1 points for a certain landmark i receives, may be regarded as aconfidence rate that the particular point is the searched for landmark.

Step 2. Minimal Cost Path Optimization

According to the first step, each landmark may separately be displacedover 2n_(s)+1 possible locations. The confidence rates N_(i,j) can beput in a matrix (218) as follows $R = {\begin{bmatrix}N_{1,1} & \cdots & N_{1,{{2n_{s}} + 1}} \\\vdots & N_{i,j} & \vdots \\N_{n,1} & \cdots & N_{n,{{2n_{s}} + 1}}\end{bmatrix}\quad{with}}$${\sum\limits_{j = 1}^{{2n_{s}} + 1}N_{i,j}} = N$ ∀i = 1, …  , n.

Updating the contour involves choosing one element in each row of thematrix C; the concatenation of all elements forming a path through thematrix, such that a cost criterion is optimized (i.e. minimized), is aminimal cost path (MCP) search.

A logical choice would be to take the element with highest confidence ineach row; this approach has the drawback that it does not take intoaccount the presence of outliers. Therefore, it may include abruptchanges in the path. This drawback is overcome by imposing a smoothnessconstraint that excludes the abrupt changes.

To this purpose, first a cost matrix is constructed as follows:${C = \begin{bmatrix}{1/N_{1,1}^{m}} & \cdots & {1/N_{1,{{2n_{s}} + 1}}} \\\vdots & {1/N_{i,j}^{m}} & \vdots \\{1/N_{n,1}^{m}} & \cdots & {1/N_{n,{{2n_{s}} + 1}}^{m}}\end{bmatrix}},$

with the power m in the denominator being a positive number thatinfluences the computed path (a larger m yielding smaller costcontribution for a point and hence steering the path towards thispoint). Second, a cost function is constructed as follows:${{J\left( {j_{1},j_{2},\ldots\quad,j_{n}} \right)} = {{\sum\limits_{i = 1}^{n}C_{i,j_{i}}} + {\alpha{\sum\limits_{i = 1}^{n}{{\delta_{i,j_{i}} - \delta_{{i - 1},j_{i - 1}}}}}}}},$

with δ_(i,j) ^(i) the displacement of point i towards its target pointj_(i) along the normal to the shape and α a weight factor. The secondterm in this equation is the smoothness constraint that penalizes targetpoint configurations for which displacements of one target point w.r.t.the displacement of the previous target point vary abruptly along thecontour. The cost function J is optimized (220), using e.g. a techniqueknown in the prior art as dynamic programming, described in e.g. D. H.Ballard, C. M. Brown, Computer Vision, Englewood Cliffs, Prentice HallInc. 1982, pp 137-143, resulting in the target points (j^(*) ₁ , . . . ,j^(*) _(n)), that minimize J as a minimal cost path (MCP).

Alternatively to the matrix data structure, a layered graph may be usedalso. Each layer represents the array of perpendicularly sampled imagepoints at a landmark. The graph is augmented with a set of arcs from thelast layer back to the first. The dynamic programming technique extractsthe boundary as an optimal closed path through the graph.

Step 3. Unweighted and Weighted PCA Contour Fitting (223)

In the previous step, each landmark i is given a new location j^(*)_(i). The sequence (j^(*) ₁ , . . . j^(*) _(n)) represents a contourthat may be fitted in the shape model x≈ x+Φb by substituting it for xand deriving b by solving Φb=x− x. Fitting x in the shape model meansfinding a b such that x+Φb approximates x with small errors in the highconfidence landmarks.

Un-weighted PCA Contour Fitting

Since Φ∈

there are more equations than unknowns, and the system isover-determined, having no exact solution. Therefore, ∥Φb−(x− x)∥_(p)must be minimized for some suitable choice of p. Different norms renderdifferent optimum solutions. In particular the case of p=2, the leastsquares problem, is more tractable, and is solved in the prior art bymethods based on normal equations and the QR factorization. Denoting thedifference between actual point sequence and the fitted point sequenceas the fit error e=x−( x+Φb), the problem is to seek the b to minimizee^(T)e. Substituting for e, this optimization problem becomes${\min\limits_{b}\quad{J(b)}} = {{\min\limits_{b}\quad{{\mathbb{e}}^{T}e}} = {\min\limits_{b}{\left( {\left( {x - \overset{\_}{x}} \right)^{T} - {b^{T}\Phi^{T}}} \right){\left( {\left( {x - \overset{\_}{x}} \right) - {b\quad\Phi}} \right).}}}}$

J(b)is minimized if ∇J(b)=0. The gradient ∇J(b) of J(b) can be written∇J(b)=−2Φ^(T)W²(x− x)+2Φ^(T)W² Φb.

Hence the optimization problem leads to the equations

∇J(b)=0

Φ^(T)Φb=Φ^(T)(x− x,

known as the normal equations. These are solved for b asb=(Φ^(T)Φ)⁻¹Φ^(T)(x− x)

This equation projects the data vector x of the image space into the PCAmodel space and returns the approximation coefficients b. Theapproximation result b may be projected back in image space byevaluating the shape model x= x+Φb.

Weighted PCA Contour Fitting (222)

The corresponding number of votes N_(i,j) ^(*) ^(i) may be used asweights when fitting the contour into the shape model. As a consequence,the PCA model will tend to give a high priority to the high confidencepoints. Points with a low confidence rate will lessen their influence onthe contour that the model generates. The optimization problem nowincludes minimizing the weighted quadratic error (We)^(T)(We). Using asimilar approach with normal equations as in the un-weighted case, thisproblem leads to the solution for b expressed asb=(Φ^(T)W²Φ)⁻¹Φ^(T)W²(x− x).

This equation is evaluated to project x in the PCA space and may beprojected back using the shape model x= x+Φb.

Steps 1 to 3 are iterated until a stop criterion is fulfilled, e.g. whenthe positions of the landmarks do not change substantially betweensuccessive iterations.

Both the un-weighted and weighted PCA fitting yields the finalsegmentation result (224) as a set of landmark points x in the image. Acontinuous spline curve may be interpolated through these points torepresent the segmentation analytically.

Modeling System II (FIG. 3)

Gray Level Appearance Model

The gray level appearance model (118′) captures the spatial intensitydistribution in the neighborhood of landmarks. In this embodiment, thespatial model departs from circular profiles sampled in the neighborhoodof a landmark. The number n of landmarks that is used to accuratelydescribe a lung field outline may be as low as 14, e.g. landmark 1denotes the uppermost point of each lung field, landmark 7 is positionedat the cross section of either the left or right heart shadow with theirrespective hemi-diaphragm, landmark 9 denoting the costophrenic angle.The other landmarks are equidistantly positioned between these majorlandmarks. Obviously, a higher number of landmarks may be used also inthe computational scheme according to the present invention.

Circular Profiles

An alternative way for selecting the points in the neighborhood of alandmark is to sample points on a circle with center at the landmark andat least one radius (312). An example of a circular sampling is given inFIG. 9. If the landmark is located at (x₀,y₀), the gray value functionof the image is sampled at radius r_(c) at the points$\left( {{x_{0} + {r_{c}\cos\frac{2\pi}{n_{c}}\left( {i - 1} \right)}},{y_{0} + {r_{c}\sin\frac{2\pi}{n_{c}}\left( {i - 1} \right)}}} \right)$i = 1, …  , n_(c).

The samples are put in a profile vector of length n_(c). Suitablechoices of n_(c) are 12, 8, 6, 4 and 3 (corresponding to 30, 45, 60, 90and 120 degree angular spacing). Multiple circular sub-samplings at aset of radii may be considered simultaneously. Suitable radii r_(c),expressed dimensionless as a fraction of the image size, are 0.05 and0.025.

This scheme has the advantage over the linear profiles that it capturesthe 2D structure around the neighborhood. Specific anatomical landmarkssuch as the costophrenic corner point or the intersection of the leftheart shadow with the left diaphragm are better modeled by circularprofiles than linear profiles for they have discontinuous firstderivative of the lung field contour. Landmarks along the lung fieldcontour in the vicinity of these specific anatomical landmarks may alsoambiguously be represented by a linear profile if its length is toolong, for then the profile may cross the juxtaposed lung field border.

The circular profile may be normalized such that the sum of absolutevalues of the elements in the profile is 1. The degree of sub-samplingmay range from no sub-sampling (yielding a re-sampling at originalresolution of the image) to very sub-sampled, e.g. only 4 points on themain axes in the case of a circular profile.

Gray-Level Feature Images

The image structure around a landmark is captured in a number N ofmathematical features as outlined in the summary of the invention. Thefeature images (314) are obtained by computing a number of moments ofthe local histogram of derivative images in a window with width α aroundeach location x₀. In a preferred embodiment the first and secondhistogram moments, all derivatives up to second-order(L,L_(x),L_(y),L_(xx),L_(yy),L_(xy)) and 5 inner scales (σ=0.5,1,2,4,8pixels) are computed, amounting to a total of N=2×6×5=60 feature images.The window scale to compute the histogram is in accordance to the innerscale, i.e. α=2σ. The set of feature images can be regarded as agray-level decomposition of the gray-level function in the neighborhoodof a landmark.

Gray-Level Cost

The Mahalanobis distance between a new profile g_(i,j) sampled in a newfeature image j for a certain landmark i is computed ash_(i,j)(p)=(g_(i,j)(p)− g _(i,j))^(T)S⁻¹ _(i,j)(g_(i,j)(p)− g _(i,j))

A smaller Mahalanobis distance means a larger probability that theprofile g_(i,j)(p) originates from the Gaussian distribution with mean g_(i,j) and covariance S_(i,j). Hence, the Mahalanobis distance may beused as a gray-level cost measure, denoted as h_(i,j)(p). This costmeasure is a function of the location p in a given image. The locationthat most likely is the true location for the landmark i, in accordanceto feature j, is the point p that minimizes h_(i,j)(p). This gray-levelcost can thus be defined for each feature j.

A total gray-level cost is obtained by combining all gray-level costsfor all features j in a sum${h_{i}(p)} = {\sum\limits_{j = 1}^{N}{\left( {{g_{i,j}(p)} - {\overset{\_}{g}}_{i,j}} \right)^{T}{S_{i,j}^{- 1}\left( {{g_{i,j}(p)} - {\overset{\_}{g}}_{i,j}} \right)}}}$

reflecting the similarity between the gray-level pattern at p and theexpected gray-level pattern at landmark i. The location most likely tocoincide with landmark i, in accordance to the gray-level appearance, isthe solution of the optimization problem$p_{o} = {\arg\quad{\min\limits_{p}{{h_{i}(p)}.}}}$

To construct the gray-level appearance model (118′), the distributions (g _(i,j),S_(i,j)) of the feature profiles g_(i,j) (for all landmarks i=1, . . . , n, for all feature images j=1 , . . . , N have to be estimatedfrom the training images.

Shape Model (126′)

Whereas a gray-level cost validates a gray-level pattern or a pattern offeatures, a shape cost is defined to validate a connection between twosuccessive landmarks.

The curve outline, one for each of the lung fields, is described by npoints (landmarks). These landmarks are manually determined in a set ofs training images, yielding a sequence of points (x₁,y₁) . . .(x_(n),y_(n))=(p₁ , . . . , p_(n)). These coordinate tuples aresubsequently positioned in a vector x=(x₁,y₁ . . . , x_(n), y_(n))^(T)(320), representing the shape. Considering a pair (p_(i),p_(i+1)) ofestimated positions of successive landmarks. A shape cost is assigned tothe vector v_(i)=p_(i+1)−p_(i), reflecting the plausibility of v_(i)w.r.t. its probability distribution. The distributions of v₁,v₂ , . . ., v_(n), assumed to have normal distributions around their mean, areestimated from the training shapes. The mean (322) and covariance (324)of v_(i) are noted as v _(i) and S_(v) ^(i) respectively. An example ofthese vector distributions is given in FIG. 10. The vector distributionscan be regarded as a first-order (tangential) geometric decomposition ofthe shape represented by the set of landmarks.

A novel shape cost validating a connection vector v_(i) between twosuccessive landmarks p_(i),p_(i−1) (a vector in the plane being fullycharacterized by its orientation and its length) is the Mahalanobisdistance between v_(i) and its mean v _(i):ƒ_(i)(v_(i))=(v_(i)− v _(i))^(T)S⁻¹ _(v) ^(i) (v_(i)− v _(i)).

A connection, which is unlikely to happen, because it has largeMahalanobis distance, will get a high shape cost. On the other hand, azero cost will be assigned to a connection that equals the expectedvalue.

Segmentation System II (FIG. 4)

Knowledge about the gray-level appearance of the object in the image andabout the shape of the object to be segmented is acquired during thetraining phase. This knowledge is subsequently used to segment theobject in a new image (418) according to the following steps.

Step 1. Rectangular Search Grids Construction (410).

A rectangular search grid for each landmark i is constructed toconstrain the possible locations of the landmark. Each rectangular gridis characterized by the parameters of geometric extent and grid spacing.

The grid spacing r_(g) should be in accordance with the radius r_(c) ofthe circular profiles (312) of the gray-level appearance model. Asuitable relationship is r_(g) =ƒ_(g)r_(c) with r_(c) the radius andƒ_(g) a fixed fraction. A typical fraction for the grid spacing isƒ_(g)=0.5.

The grid extent and grid position is landmark specific and isrepresented by x_(min), x_(max), y_(min) and y_(max). The true (unknown)landmark location (x^(*),y^(*)) is deemed to lie between the lower andupper borders:x_(min)<x^(*)<x_(max)y_(min)<y^(*)<y_(max)

These conditions can only be guaranteed if the grid spans the wholeimage area (i.e. x_(min)=y_(min)=0 and x_(max)=y_(max)1). A moreefficient approach is to constrain the search to a relevant part of theimage. Assuming that the probability distribution p_(x)(x) of thex-coordinate of the landmark is normal with mean x and standarddeviation σ_(x) (estimated from the training shapes). The grid bordersx_(min) and x_(max), jointly defining the grid extent in thex-direction, are chosen such that∫_(x_(min))^(x_(max))p_(x)(x)  𝕕x = f_(c)

with ƒ_(c) some fraction representing the probability that the conditionx_(min)<x^(*)<x_(max) is true. A suitable value for the fraction ƒ_(c)is 0.995. If the requirement is imposed that the extent is symmetricaround the mean x, i.e. x_(max)− x= x−x_(min), the upper border x_(max)is the value that satisfies${\frac{1}{\sigma_{x}\sqrt{2\pi}}{\int_{\overset{\_}{x}}^{x_{\max}}{{\exp\left( {- \frac{\left( {x - \overset{\_}{x}} \right)^{2}}{2\sigma_{x}^{2}}} \right)}\quad{\mathbb{d}x}}}} = {\frac{1}{2}{f_{c}.}}$

The lower border is then x_(min)=2 x−x_(max). The borders for the ycoordinate are obtained analogously.

Step 2. Gray-Level Cost Matrix Construction

The gray-level appearance model (126′) is used to find proper locationsfor each landmark. The m best locations of landmark i, in accordancewith the gray level appearance model, are selected as follows.

-   First, a rectangular grid covering a relevant part of the image area    around the expected location of each landmark i is defined,    according to the previous step.-   Secondly, for each grid point, the total gray-level cost h_(i)(p) is    computed. The points p_(i,1),p_(i,2) , . . . p_(i,m) corresponding    to the m lowest total gray-level costs are selected.-   Thirdly, a gray-level cost matrix C (414) is constructed so that    each row i contains the costs of the m most likely locations of    landmark i: $C = {\begin{bmatrix}    {h_{1}\left( p_{1,1} \right)} & \cdots & {h_{1}\left( p_{1,k} \right)} & \cdots & {h_{1}\left( p_{1,m} \right)} \\    \vdots & \quad & \vdots & \quad & \vdots \\    {h_{\quad i}\left( p_{\quad{i,\quad 1}} \right)} & \cdots & {h_{i}\left( p_{i,k} \right)} & \cdots & {h_{i}\left( p_{i,m} \right)} \\    \vdots & \quad & \vdots & \quad & \vdots \\    {h_{n}\left( p_{n,1} \right)} & \cdots & {h_{n}\left( p_{n,k} \right)} & \cdots & {h_{n}\left( p_{n,m} \right)}    \end{bmatrix}.}$

A typical number m of best points per landmark ranges from 10 to 40.

As the m of best points per landmark are selected independently from theset of m of best points for a neighboring landmark, the situation mayarise that one or more of these points are nearer to a non-neighboringlandmark. These points will likely be neglected in the final contoursegmentation by accounting for the shape cost in the next step.

Step 3. Minimal Cost Path Construction (416)

Determining the contour that segments the object is reduced to finding apath from top to bottom through the matrix C by selecting one elementper row. Denoting the index of the selected element in row i as k_(i),the curve becomes the point sequence {p_(1,k) ¹ ,p_(2,k) ² , . . .p_(n,k) ^(n) }. The optimal path (k^(*) ₁,k^(*) ₂ , . . . , k^(*) _(n))is the path that minimizes a cost function J(k₁ , . . . , k_(n)):$\left( {k_{1}^{*},k_{2}^{*},\ldots\quad,k_{n}^{*}} \right) = {\arg\quad{\min\limits_{k_{1},\ldots\quad,k_{n}}{{J\left( {k_{1},\ldots\quad,k_{n}} \right)}.}}}$

The models introduced above admit a number of cost measures. Consideringtwo successive landmarks p_(i,k) ^(i) and p_(i+1,k) ^(i+1) , a costcomponent (412) according to the gray-level appearance model and a costcomponent according to the shape model are:

-   The gray level costs h_(i)(p_(i,k) ^(i) ) and h_(i+1)(p_(i+1,k)    ^(i+1) ) corresponding to the landmarks pink and p_(i,k) ^(i) and    p_(i+1,k) ^(i+1) ;-   The shape cost ƒ_(i)(p_(i+1,k) ^(i+1) −p_(i,k) ^(i) ), which    validates the plausibility of the connection from p_(i,k) ^(i) to    p_(i+1,k) ^(i+1) .

Both cost types can be incorporated in the cost function J(k₁ , . . .,k_(n)) by means of a weighted sum of an overall gray-level cost and anoverall shape cost. The overall gray-level cost with weight factor γ₁ isthe sum of the landmark individual gray-level costs h_(i)(p_(i,k)_(i)),i=1 , . . . , n. The overall shape cost with weight factor γ₂ isthe sum of the shape costs ƒ_(i)(p_(i−1,k) ^(i+1) −p_(i,k) ^(i) ),i=1 ,. . . , n. Hence the cost function J(k₁ , . . . , k_(n)) becomes${J\left( {k_{1},\ldots\quad,k_{n}} \right)} = {{\gamma_{1}{\sum\limits_{i = 1}^{n}\quad{h_{i}\left( p_{i,k_{i}} \right)}}} + {\gamma_{2}{\sum\limits_{i = 1}^{n - 1}\quad{f_{i}\left( {p_{{i + 1},k_{i + 1}} - p_{i,k_{i}}} \right)}}} + {\gamma_{2}{{f_{n}\left( {p_{1,k_{1}} - p_{n,k_{n}}} \right)}.}}}$

The optimal path (k^(*) ₁,k^(*) ₂ . . . , k^(*) _(n)) that minimizesJ(k₁ , . . . , k_(n)) is called a Minimal Cost Path. The optimal path(416) is computed using dynamic programming techniques known in theprior art, such as introduced in G. Behiels et al., Evaluation of imagefeatures and search strategies for segmentation of bone structures usingactive shape models, Medical Image Analysis, 6(1):47-62, 2002. A typicalweight factor for the gray-level cost is γ₁=1, and for the shape costγ₂=0.25.

FIG. 13 shows the automatic segmentation of three thorax imagesaccording to the system outlined before, and compares it to the manualsegmentation performed by an experienced radiologist. The 14 auto-placedlandmarks on each lung field are shown, along with a continuous splineinterpolation through them.

Computational Speed-up Refinements

The segmentation method including steps 1-3 can be applied to segmentonly part of the outline of one or more anatomic entities. To thispurpose the rectangular search grid is only constructed for a subset ofconsecutive landmarks. The gray level cost matrix will comprise a numberof rows equal to the number of retained landmarks. The minimum cost pathconstruction will minimize a cost function that only comprises a graylevel term and a shape term for the retained landmarks. The resultingsegmentation traces only that outline part that is covered by theretained landmarks. The obvious advantage of partial segmentation iscomputational speed when one is only interested in part of the anatomicentity outline, or when only some specific landmarks need to be locatedto obtain specific measurements.

Another algorithmic refinement to speed up the computation of thesegmentation is to use a coarse-to-fine multi-resolution approach. Thegray-level appearance and shape models are constructed on amulti-resolution representation for each landmark, introducing thespatial scale as a new dimension. Two alternatives are implemented tothis purpose, both applying the segmentation concept over successivelyfiner resolution levels.

In one alternative, during segmentation, the number of points in aprofile, the number of points in the search grid and the number oflandmark points along the contour are decreased in the original image.The positional precision of the landmarks is then increased bysuccessively increasing the resolution of the search grid in thevicinity of the previous solution (the total number of considered pointsin the profile and search grid may thus be held constant at eachiteration), and applying the finer resolution gray-level and shapemodels.

In another alternative, the multi-resolution approach may be implementedusing a number of levels of a Gaussian and Laplacian pyramid and theirassociated derivative feature images. The initial position of thelandmarks is determined on a coarse level, and is tracked to its finalposition on the finest level using the intermediate positions of theprevious level as the starting position for the next level. Because thecoarse resolution images contain fewer details, the search space willlikely contain less false minima which will enable faster optimizationat a given level and faster locating the landmarks towards their finalposition.

Training the Segmentation Models

A number of thoracic images is collected from an image repository, andpresented sequentially in a Graphical User Interface. The followingsteps are performed to build the segmentation models:

-   An experienced user manually segments the lung fields in each    displayed image, by indicating—using left mouse clicks—a number of    points along the contour that delineates the lung field. These    points need not be spaced equidistantly along the lung field    outline; the only requirement is that the point set collectively    approximates the lung field to a sufficiently high degree. To assess    the anatomic fit, a spline curve is continually interpolated through    the manually positioned points determined so far, until the curve is    closed from the last point to the first point by a right mouse    click. Manual segmentation adjustments can be made by dragging an    individual point towards a new position. The resulting outline is    again assessed with respect to anatomical correctness.-   Next, a number of landmarks will be auto-placed on the manual    segmented lung field contour. In order to achieve that identical    landmarks on all images of the training set map on each other, the    user positions a few number of easily discernible landmarks, the    other landmarks are obtained by equidistantly placing points on the    lung field contour. In the case of lung field outlines, a number of    easily discernable landmarks are the topmost lung field point, the    point denoting the costophrenic angle and the junction of heart    shadow and hemi-diaphragm, a sub-total of three. Next a number of    total landmarks are chosen, a suitable choice ranging from 14 to 40    for example. In the 14 landmark case, the points p₁, p₇ and p₉    represent the three fixed landmarks, five points are evenly    distributed between p₁, and p₇, one point between p₇ and p₉, and    another five points between p₉ and p₁. This step is performed    separately on left and right lung fields.-   Next, parameters pertaining to the training phase are asked: (a)    image size for training (and segmentation), a typical value is 256    or 512; (b) fraction of the shape variance to be explained by the    Principal Components Analysis (PCA), a typical value being    0.9999; (c) number of points in the profile, either linear or    circular, a typical value for the circular profile being 12, 8, 6, 4    or 3; (d) radius of the circular profile as a fraction of the image    dimension, typical values being 0.04, 0.03, 0.02.-   Training of the gray-level appearance model: (a) decomposition of    the gray-level function around the landmarks, i.e. computation of N    feature images e.g. the local histogram moments of gray-value    derivatives as outlined before, and (b) collection of the gray value    distribution of the linear or circular profiles at the landmarks    p_(i), i.e. computation of g _(i,j),i=1 . . . n,j=1 . . . N and the    covariance matrices S_(i,j),i=1 . . . n,j=1 . . . N (i.e. for each    landmark i in feature image j). This step generates the statistical    gray-level appearance knowledge.-   Training of the shape model: (a) computation of geometric    decomposition of the zeroth order (positional) information contained    in the landmarks, i.e. computation of the mean shape x and the t    principal modes of variation (eigenshapes) arranged columnwise in a    matrix b; (b) computation of the vector distributions (first order    tangential information) contained in the connection sequence of    landmarks, i.e. computation of the mean v _(i) and covariance S_(v)    ^(i) of v_(i). This step generates the statistical shape knowledge.-   Storage of the statistical gray-level appearance and shape    knowledge, e.g. to be used to segment the lung fields in a new image    according to the model-based segmentation sub-systems given above.    Application of the Model-based Segmentation in 2D Images to the    Computation of the Cardio-Thoracic Ratio (CTR)

An application of the automated lung field segmentation is thecomputation of the Cardiothoracic Ratio (CTR). The CTR (FIG. 6) isdefined as the ratio of the transverse diameter of the heart to theinternal diameter of the thorax (ID): ${CTR} = \frac{{MLD} + {MRD}}{ID}$with MLD the maximum transverse diameter on the left side of the heartand MRD the maximum transverse diameter on the right side of the heart.This index is an important clinical parameter, which varies for an adultbetween 39% and 50% with an average of about 45%. A cardiothoracic indexhigher than 50% is considered abnormal. Possible causes are cardiacfailure, pericardial effusion and left or right ventricular hypertrophy.It is possible to compute the cardiothoracic ratio automatically, usingthe automatic lung field segmentation as disclosed in the presentinvention.

Referring to FIG. 8, the characteristic point defining MRD is obtainedby selecting the Cartesian leftmost point on the fitted contour segmentbetween landmarks p1 and p9 of the right lung segmentation; thecharacteristic point defining MLD is obtained by selecting the Cartesianrightmost point on the fitted contour segment between landmarks p1 andp7 of the left lung segmentation. The sum MLD+MRD is obtained bysubtracting the column coordinates of these characteristic points andtaking the absolute value. Similarly, the ID is obtained by selectingthe Cartesian leftmost point on the fitted contour segment betweenlandmarks p1 and p7 of the right lung segmentation and the Cartesianrightmost point on the fitted contour segment between landmarks p1 andp9 of the left lung segmentation, subtracting the column coordinates ofthese characteristic points and taking the absolute value.

FIG. 14 shows the computation of the characteristic points on thesegmented lung fields, and the computation of the cardiothoracic ratiousing (a) the manual lung field segmentations resulting in aCTR_(man)=0.47; (b) using the automatically derived lung fieldsegmentations according to the present invention resulting in aCTR_(automatic)=0.45.

Application of the Model-based Segmentation and Measurement System toOther Body Parts

The spatial extent of the search grid for each of the landmarks in thelung field segmentation is derived on the basis of the positions of allsimilar landmarks in a set of training thorax images. The concept ofusing a search grid for constraining the candidate positions for a givenlandmark can be extended for other body parts with that may have widerpositional, rotational and size variation in the image than that of alung field in a thorax image, and that do occupy the full image area asopposed to lung fields in a thorax image.

To allow for such positional, rotational and size variation of the bodypart, the concept of anchor point mapping of the search grids usingmethods disclosed in European patent application 04076454, entitled“Method for automatically mapping of geometric objects in digitalmedical images”, may be applied in conjunction with the presentinvention. The refined method then becomes particular interesting tosegment bony body parts in orthopedic radiographs because they arecharacterized by a relatively constant shape in the image, and they canbe anchored to well manifested bony landmarks in the image.

For example, in a pelvis, hip or full leg examination, the femoraloutlines can be anchored to well known landmarks such as the tip of thegreater trochanter, the knee center and the femoral head center. Theseanchor points are used to establish an affine transformation betweenmodel anchor points and the associated anchor points selected in theactual image. Because a search grid has a collection of candidate pointsarranged in a rectangular lattice around a landmark point on thesegmentation model contour in the model image, each of the constituentcandidate points is mapped in turn in the actual image by applying thetransformation to the model points' coordinates. In this way, the searchgrid for a given landmark point is reconstructed around the most likelyposition of the point. The optimization process then proceeds in themanner disclosed above by optimizing a combined gray value and shapecost for a path through a selected combination of candidate locations,one on each search grid. The path with most optimal cost is the finalsegmentation of the body part. Examples of other bones are hand andupper extremity bones, foot and lower extremity bones, pelvis and spinalvertebrae and structures.

Other body parts that are amenable to this type of landmark-basedsegmentation are soft-tissue organs in 3D CT of MR images. On the onehand, a slice-by-slice based approach may be taken here, that determinesthe segmented contour of the anatomy on each slice, and that combinesthe result of all slices into a 3D segmented surface. On the other hand,a fully 3D approach, extending the current concept on a 3D volume, maybe adopted, that builds and applies models for 3D landmarks confinedwithin 3D search grids. Kidneys, lungs, heart, liver, stomach, spleen,prostate and brain structures are examples of organs to which the methodof the present invention is applicable.

In all application examples mentioned, measurement points may be basedon the position of the landmarks on the segmentation or on a combinationof several landmarks (such as a midpoint of two juxtaposed points oneither side of the cortex of an elongated bone, a pair of whichdetermines the bone's anatomical axis), and the landmark points may besubsequently fed into a measurement system such as disclosed in EP A1349098.

1. A method of segmenting an anatomic entity in a digital medical imagecomprising the following steps: defining, for each of a number oflandmarks in said image, an initial position of said landmark, samplinga neighborhood around said initial position, said neighborhoodcomprising a number of candidate locations of said landmark, associatinga cost with each of said candidate locations, optimizing a cost functionexpressing a weighted sum of overall gray level cost and overall shapecost for all candidate locations, defining a segmented anatomic entityas a path through a selected combination of said candidate locations forwhich combination said cost function is optimized, wherein said costassociated with a candidate location is the total gray level costobtained by summing gray level costs of said candidate location in everyfeature image of a set of feature images calculated for said medicalimage, said gray level costs being expressed as the Mahalanobis distancedefined as the distance between the gray value profile at said candidatelocation with a corresponding mean profile weighted by the inverse of acorresponding covariance matrix, said mean profile and covariance matrixretrieved from a gray value model associated with said anatomic entity.2. A method according to claim 1 wherein said gray value model isobtained by the steps of sampling a manually segmented outline of saidanatomic entity at a number of landmark points; sampling a number ofpoints in a neighborhood of each of said landmark points; for eachlandmark point arranging sampled points in a neighborhood around alandmark point in a profile; computing at least one feature image;computing a mean profile for each landmark point and for each featureimage; computing a covariance matrix of the profiles for each landmarkpoint and each feature image; identifying said mean profile andcovariance matrix as gray value model of said anatomic entity.
 3. Amethod according to 1 wherein a feature image is a derivative image at apredefined scale σ or an image representation comprising mathematicalmoments of the local histogram of said image or derivative image in awindow with width α around the location of a landmark.
 4. A methodaccording to claim 1 wherein said overall gray level cost is acombination of all individual gray level costs of a profile in a numberof feature images whereby said individual cost is the Mahalanobisdistance defined as the distance between said profile and the meanprofile of that feature weighted with the inverse of the covariancematrix, said mean profile and covariance matrix being retrieved from agray value model of said anatomic entity.
 5. A method according to claim4 wherein said gray value model is obtained by the steps of sampling amanually segmented outline of said anatomic entity at a number oflandmark points; sampling a number of points in a neighborhood of eachof said landmark points; for each landmark point arranging sampledpoints in a neighborhood around a landmark point in a profile; computingat least one feature image; computing a mean profile for each landmarkpoint and for each feature image; computing a covariance matrix of theprofiles for each landmark point and each feature image; identifyingsaid mean profile and covariance matrix as gray value model of saidanatomic entity.
 6. A method according to claim 1 wherein said overallshape cost is a combination of all individual costs of connectionvectors, said connection vectors connecting successive landmarks, forall landmarks whereby said individual shape cost is expressed by theMahalanobis distance defined as the distance of a connection vectorbetween two successive landmarks to the mean connection vector weightedwith the inverse of the covariance matrix, said mean connection vectorand covariance matrix being retrieved from a shape model of saidanatomic entity.
 7. A method according to claim 6 wherein said shapemodel is obtained by the steps of sampling a manually segmented outlineof said anatomic entity at a number of landmark points; computingconnection vectors or connection vector differences between successivelandmark points; computing a mean connection vector or mean connectionvector difference for successive pairs of landmark points; computing acovariance matrix of connection vectors or connection vectordifferences; identifying said mean connection vector and covariancematrix as a geometric model of said anatomic entity.
 8. A methodaccording to claim 1 wherein said cost function is optimized by dynamicprogramming.
 9. A method according to claim 1 wherein said number ofcandidate locations is reduced by selecting those with minimal totalgray value cost, said total gray value cost being the sum over allfeature images of all gray value costs for the candidate position insaid neighborhood, said gray value cost expressed as a Mahalanobisdistance defined as the distance of a gray value profile at thecandidate location to a corresponding mean profile weighted with theinverse of a covariance matrix, said mean profile and covariance matrixretrieved from a gray value model associated with said anatomic entity.10. A method according to claim 1 wherein said neighborhood comprises arectangular grid of sampling points.
 11. A method according to claim 1wherein said neighbourhood comprises a circular profile.
 12. A methodaccording to claim 1 wherein said gray value model and said geometricmodel are constructed from and applied to a multi-resolutionrepresentation of the image.
 13. A method of measuring an anatomicentity in a digital medical image comprising the steps of calculating asegmentation of said anatomic entity using a method according to claim1, calculating or selecting characteristic points based on thesegmentation of said anatomic entity, calculating measurement objectsbased on said characteristic points, deriving measurements of saidanatomic entity on the basis of said measurement objects.
 14. A methodof measuring an anatomic entity in a digital medical image comprisingthe steps of calculating a segmentation of said anatomic entity using amethod according to claim 1, identifying said landmarks as measurementpoints, calculating measurement objects based on said measurementpoints, deriving measurements of said anatomic entity on the basis ofsaid measurement objects.
 15. A computer program product embodied in acomputer readable medium for performing the steps of the method ofclaim
 1. 16. A computer readable medium comprising computer executableprogram code for performing the steps of the method of claim
 1. 17. Acomputer program product embodied in a computer readable medium forperforming the steps of the method of claim
 13. 18. A computer readablemedium comprising computer executable program code for performing thesteps of the method of claim
 13. 19. A computer program product embodiedin a computer readable medium for performing the steps of the method ofclaim
 14. 20. A computer readable medium comprising computer executableprogram code for performing the steps of the method of claim
 14. 21. Amethod of segmenting an anatomic entity in a digital medical imagecomprising the following steps: defining, for each of a number oflandmarks in said image, an initial position of said landmark, samplinga neighborhood around said initial position, said neighborhoodcomprising a number of candidate locations of said landmark, associatinga cost with each of said candidate locations, optimizing a cost functionexpressing a weighted sum of overall gray level cost and overall shapecost for all candidate locations, defining a segmented anatomic entityas a path through a selected combination of said candidate locations forwhich combination said cost function is optimized, wherein said costassociated with a candidate location is the total gray level costobtained by summing gray level costs of said candidate location in everyfeature image of a set of feature images calculated for said medicalimage, said gray level costs being expressed as the Mahalanobis distancedefined as the distance between the gray value profile at said candidatelocation with a corresponding mean profile weighted by the inverse of acorresponding covariance matrix, said mean profile and covariance matrixretrieved from a gray value model associated with said anatomic entity,and wherein said neighborhood comprises a circular profile.
 22. A methodof segmenting an anatomic entity in a digital medical image comprisingthe following steps: defining, for each of a number of landmarks in saidimage, an initial position of said landmark, sampling a neighborhoodaround said initial position, said neighborhood comprising a number ofcandidate locations of said landmark, associating a cost with each ofsaid candidate locations, optimizing a cost function expressing aweighted sum of overall gray level cost and overall shape cost for allcandidate locations, defining a segmented anatomic entity as a paththrough a selected combination of said candidate locations for whichcombination said cost function is optimized, and wherein said overallshape cost is a combination of all individual costs of connectionvectors, said connection vectors connecting successive landmarks, forall landmarks whereby said individual shape cost is expressed by theMahalanobis distance defined as the distance of a connection vectorbetween two successive landmarks to the mean connection vector weightedwith the inverse of the covariance matrix, said mean connection vectorand covariance matrix being retrieved from a shape model of saidanatomic entity.